%% replicate_fig6a.m
clear; clc; close all;

%% 参数设置
a     = 0.6;      
b     = 0.5;
alpha = 2.1;      
beta  = 2;
gamma = 2;
Im    = 0.15;
F     = 0.15;

%% 积分时间设置
T_final = 3500;  % 总积分时间
T_trans = 3000;  % 舍去瞬态
tspan   = [0, T_final];

% 初始条件
Y0 = [0; 0; 0];  % [q; sigma; phi]

% ODE45 求解选项
opts = odeset('RelTol',1e-9,'AbsTol',1e-9);

%% 数值积分
[t, Y] = ode45(@(t,Y) cnn(t, Y, a, b, alpha, beta, gamma, Im, F), tspan, Y0, opts);

% 取出各状态
q     = Y(:,1);
sigma = Y(:,2);
phi   = Y(:,3);

% 计算 Vx = (alpha + beta*sigma).*q
Vx = (alpha + beta*sigma).*q;

%% 舍去瞬态，保留 [T_trans, T_final] 部分
idx     = (t >= T_trans);
t_plot  = t(idx);
phi_plot= phi(idx);
Vx_plot = Vx(idx);

%% 绘图
figure;

% 子图1：相平面图，占图像左边 35% 的宽度
subplot(1, 2, 1);
plot(phi_plot, Vx_plot, 'r', 'LineWidth', 1.5);
xlabel('\phi','FontSize',12);
ylabel('V_x','FontSize',12);
title('相平面图：\phi vs. V_x','FontSize',12);
grid on;

% 子图2：时间波形图，占图像右侧 50% 的宽度
subplot(1, 2, 2);
plot(t_plot, Vx_plot, 'b', 'LineWidth', 1.5);
xlabel('时间 t','FontSize',12);
ylabel('V_x','FontSize',12);
title('V_x 随时间变化波形','FontSize',12);
grid on;

function dYdt = cnn(t, Y, a, b, alpha, beta, gamma, Im, F)
    % 无微扰 CNN 单元模型 (5)
    % Y = [q; sigma; phi]

    q     = Y(1);
    sigma = Y(2);
    phi   = Y(3);

    % 计算 V_x = (alpha + beta*sigma) * q
    Vx = (alpha + beta*sigma)*q;

    % 外激励信号
    Iext = Im * sin(2*pi*F*t);

    % memristor 的非线性函数 G_M = phi^2 - 0.5
    G_M = phi^2 - 0.5;

    % 模型 (5) 的微分方程
    dq_dt     = Iext - G_M*Vx + a*abs(Vx) + b*Vx;
    dsigma_dt = q - sigma;
    dphi_dt   = gamma*Vx + tanh(phi) - phi;

    dYdt = [dq_dt; dsigma_dt; dphi_dt];
end
